%% fig8a_main.m
% 在 (m21, m12) 参数平面上计算“谱熵”分布，并绘制彩色热力图
% 以复刻类似论文 Fig.8(a) 的效果。

clear; clc; close all;

%% 1. 设置固定的系统参数
params.a1 = 1.5;  params.b1 = 3;  params.c1 = 1;  params.d1 = 5;  params.i1 = 0.45;
params.a2 = 0.77; params.b2 = 0.3; params.c2 = 0.8; params.i2 = 0;
params.a3 = 0.9;  params.b3 = 3;  params.c3 = 1;  params.d3 = 5;  params.i3 = 0.5;

params.m23 = 0.785;  % 其余耦合系数
params.m32 = 1.01;

eps_val = 4.87;

%% 2. 定义参数扫描范围 (m21, m12)
m21_vals = linspace(2.5, 3.5, 100);   % 示例: 2.5~3.5
m12_vals = linspace(4,   5,   100);   % 示例: 4~5

%% 3. 设置积分及谱熵计算参数
T_trans = 200;    % 去除瞬态时间
T_data  = 1000;    % 采样数据时间 (用于计算频谱)
RelTol  = 1e-6;
AbsTol  = 1e-8;

% 初始条件 (与论文一致)
X0 = [-2, 0, 0, 0, 0, 0.1, 0.3];

% 采样率 (Hz), 仅用于做谱分析时
% 若系统为无量纲，可假设 dt=1 对应 1 step, 再给定 Fs 仅用于做离散频谱
Fs = 50;   % 示例：每 1 时间单位采样 100 次 (可根据实际需求调整)

%% 4. 预分配存储谱熵
nm21 = length(m21_vals);
nm12 = length(m12_vals);
Entropy_map = zeros(nm12, nm21);

%% 5. 并行池
if isempty(gcp('nocreate'))
    parpool;
end

%% 6. 并行扫描 (m21, m12)
parfor i = 1:nm21
    m21_cur = m21_vals(i);
    rowData = zeros(1, nm12);
    
    for j = 1:nm12
        m12_cur = m12_vals(j);
        
        % 在 parfor 中复制 params
        localParams = params;
        localParams.m21 = m21_cur;
        localParams.m12 = m12_cur;
        
        % 计算谱熵
        entVal = Entropy(eps_val, localParams, X0, T_trans, T_data, Fs, RelTol, AbsTol);
        rowData(j) = entVal;
    end
    
    Entropy_map(:, i) = rowData;
end

%% 7. 绘图
figure;
imagesc(m21_vals, m12_vals, Entropy_map);
axis xy;
colorbar;
colormap jet;
xlabel('m_{21}','FontSize',12);
ylabel('m_{12}','FontSize',12);
title('Fig.8(a): Spectral Entropy in (m_{21}, m_{12}) plane','FontSize',14);
set(gca,'FontSize',12);

% 如果想加 contour 线
% hold on;
% [C,h] = contour(m21_vals, m12_vals, Entropy_map, 10, 'k');
% clabel(C,h);


function entVal = Entropy(eps_val, params, X0, T_trans, T_data, Fs, RelTol, AbsTol)
% 计算在给定 (eps_val, params) 下的 x1(t) 信号“谱熵”
%
% 输入:
%   eps_val   : FN 神经元参数 ε
%   params    : 系统参数结构体
%   X0        : 初始条件
%   T_trans   : 去除瞬态时间
%   T_data    : 实际采样数据时间 (用于计算频谱)
%   Fs        : 采样率(仅用于做谱分析)
%   RelTol, AbsTol : ode45 容忍度
%
% 输出:
%   entVal    : 计算得到的谱熵 (单位为比特或 nat 取决于对数底，下面用 log2)

    % 1) 去除瞬态
    odefun = @(t, X) hrfnhr(t, X, eps_val, params);
    options = odeset('RelTol',RelTol,'AbsTol',AbsTol);
    [~, X_sol_trans] = ode45(odefun, [0 T_trans], X0, options);
    X_start = X_sol_trans(end,:)';  % 瞬态结束状态

    % 2) 采样真实数据
    %    先定义一个离散时间序列, 每步 dt=1/Fs, 共采样 T_data
    dt = 1/Fs;
    tspan_data = 0:dt:T_data;
    [~, X_sol_data] = ode45(odefun, tspan_data, X_start, options);
    
    % 3) 提取 x1(t) 信号
    x1_data = X_sol_data(:,1);

    % 4) 计算 x1_data 的功率谱 (可用 pwelch 或 periodogram 等)
    %    这里演示使用 pwelch
    %    pwelch 返回功率谱 Pxx, 以及对应频率向量 Fxx (单位与 Fs 对应)
    %    nfft=length(x1_data), 也可自定义
    [Pxx, Fxx] = pwelch(x1_data, [], [], [], Fs);

    % 5) 归一化功率谱
    Pxx_norm = Pxx / sum(Pxx);

    % 6) 计算谱熵: -sum( p(i) * log2(p(i)) )
    %    若 p(i)=0 则 p(i)*log2(p(i))=0, 要注意数值安全
    p_nonzero = Pxx_norm(Pxx_norm>0);
    entVal = -sum(p_nonzero .* log2(p_nonzero));

end
